* D:\E\replications\BJPS2020\resultsets\figures_2_6.do

	* was C:\Dropbox\Rehm\Hacker\variability\december2017\survival.do
	* was C:\Dropbox\Rehm\Hacker\submission\BJPS_submission\ESI_by_age.do

* cd D:\E\replications\BJPS2020\
clear all
cap log close
log using ./resultsets/figures_2_6.txt, text replace
display "$S_TIME  $S_DATE"
about
ado dir

* Figure 2: Probability of experiencing at least on 25% arc income drop within time window
	* Shown in the figure are average predicted logistic hazard rates.
	* Appendix-Figure 2: Probability of experiencing at least on 50% arc income drop within time window

set more off
use if w11101>0 & w11101<. using ./rawdata/cnef/licensed/cnef_micro.dta, clear

drop if source=="UKLHS"
	// panels not long enough

destring xwaveid, gen(pid_aul)
cap drop pid
clonevar pid=x11101ll
replace pid=pid_aul if ccode==900
replace pid=x11101LL if inlist(ccode,2,365)
sum ccode pid year

isid ccode pid year
egen id=group(source ccode pid)
isid id year

clonevar age=d11101
clonevar HHsize=d11106
rename w11101 weight

* Education
gen isced3=d11108 if inrange(d11108,1,3)

clonevar preY=i11101
clonevar postY=i11102
replace postY=i11113 if ccode==2

// UKLHS has 'income last month' only
replace preY=preY*12 if source=="UKLHS"
replace postY=postY*12 if source=="UKLHS"

cap drop _merge
merge m:1 ccode year using ./rawdata/cnef/OECD_CPI.dta
drop if _merge==2
drop _merge

gen equivalence=sqrt(HHsize)
label var equivalence "Equivalence scale: sqrt(HHsize)"

* Income concepts
gen nY=postY/cpi if postY>0
gen gY=preY/cpi if preY>0
label var nY "HH yearly net income, not equivl."
label var gY "HH yearly market income, not equivl."
gen nYe=nY/equivalence
gen gYe=gY/equivalence
label var nYe "HH yearly net income, equivl."
label var gYe "HH yearly market income, equivl."

** income changes
* top- and bottom-coding at p1 and p99
foreach v of varlist nYe gYe { 
	bys ccode year: egen PCTTOP_`v'=wpctile(`v'), p(99) weights(weight)
	bys ccode year: egen PCTBOT_`v'=wpctile(`v'), p(1) weights(weight)
	replace `v'=PCTTOP_`v' if `v'>PCTTOP_`v' & `v'<.
	replace `v'=PCTBOT_`v' if `v'<PCTBOT_`v' & `v'<.
}

* Income drops
foreach v in nYe gYe {
	xtset id year, yearly
	cap drop d`v'
	gen d`v'=d1.`v'/(abs(l1.`v'+`v')/2) if `v'<. & l1.`v'<. // arc changes
	label var d`v' "Arc change of `v'"
	gen ESI25_`v'=(d`v'<=-0.25) if d`v'<.
	gen ESI50_`v'=(d`v'<=-0.5) if d`v'<.
	gen ESI10_`v'=(d`v'<=-0.1) if d`v'<.
	if "`v'"=="nYe" local foo "HH disposable income"
	else if "`v'"=="gYe" local foo "HH market income"
	foreach j in 10 25 50 {
		label var ESI`j'_`v' "`j'+% arc drop, `foo'"
	}
}

compress
save C:\Users\rehm.16\Desktop\tmp\cnef\cnef_tmp5_figure2.dta, replace
* use C:\Users\rehm.16\Desktop\tmp\cnef\cnef_tmp5_figure2.dta, clear

isid id year

* permanent income
bys ccode pid (year): egen mean_nYe=mean(nYe)
bys ccode: egen perm3_nYe=xtile(mean_nYe), nq(3)


*******************
* Figure 2
*******************
	* Hazard rates
	// https://www.iser.essex.ac.uk/resources/survival-analysis-with-stata
	// https://www.iser.essex.ac.uk/files/teaching/stephenj/ec968/pdfs/ec968st6.pdf
	
	set matsize 1000
	foreach v in ESI25_nYe ESI50_nYe ESI10_nYe { 
		cap drop min_year
		bys ccode id: egen min_year=min(year) if `v'<.
		gen n_`v'=year-min_year+1
		bys ccode id (min_year): replace min_year=min_year[1]
		cap drop sumi_`v'
		bys ccode id (year): gen sumi_`v'=sum(mi(`v')) if year>=min_year & min_year<.
		cap drop sum_`v'
		bys ccode id (year): gen sum_`v'=sum(`v')
		replace sum_`v'=. if sum_`v'>1
		replace sum_`v'=. if sumi_`v'>1
		bys ccode id (year): gen test=sum(sum_`v')
		replace sum_`v'=. if test>1
		drop test
		
		cap drop time
		gen time=n_`v'

		quietly eststo e_`v': logit sum_`v' i.time##i.ccode##i.perm
		cap drop p_`v'
		predict p_`v', p
		bysort ccode id (time): gen s_`v' = exp(sum(ln(1-p_`v'))) if p_`v'<.
		gen f_perm_`v'=1-s_`v' if p_`v'<.
		di "*** `v' ***"
		table time perm ccode [aw=weight], cont(mean f_perm_`v') f(%4.3f)
	}

	preserve
		rename perm3_nYe perm
		collapse (mean) f_perm_* [aw=weight], by(time perm ccode)
		drop if perm==.
		
		label define ccode ///
			2   "USA (PSID)" ///
			200 "Great Britain (BHPS)" ///
			225 "Switzerland (SHP)" ///
			255 "Germany (SOEP)" ///
			900 "Australia (HILDA)", modify                 
		label val ccode ccode
		
		rename f_perm_ESI* p*
		rename p??_* p??
		
		reshape long p, i(ccode perm time) j(dropsize)
		reshape wide p, i(ccode time dropsize) j(perm)
		
		label var p1 "Permanent income: low"
		label var p2 "Permanent income: medium"
		label var p3 "Permanent income: high"

		label var time "Number of years"

		note: logit sum_X i.time##i.ccode##i.perm
		note: predict p, p
		note: bysort ccode id (time): gen s = exp(sum(ln(1-p))) if p<.
		note: gen f_perm=1-s if p<.
		note: ./resultsets/figure_2_6.do
		note: Created on `= c(current_date)'
		saveold ./resultsets/figure_2.dta, replace
		
		foreach j in 10 25 50 {
			tw (connected p1 p2 p3 time if dropsize==`j' & time<=15, sort ms(oh smx none) lp(solid solid dash) msize(medlarge)) ///
				, by(ccode, r(1) colf note("")) scheme(rbn1mono) ///
				xlabel(1(2)15, angle(0)) ///
				legend(r(1) size(vsmall) region(lstyle(none))) ///
				ysize(3) xsize(5) ///
				ylabel(0(.1)1) yline(0.5) ///
				ytitle("Mean prob. of experiencing at least one `j'%+ DI income loss")

			graph save   ./resultsets/figures/figure2_`j'.gph, replace
			*graph export ./resultsets/figures/figure2_`j'.emf, replace
			graph export ./resultsets/figures/figure2_`j'.pdf, replace
		}	
		
	restore


**********************************
* Figure 6
**********************************
	foreach j in 10 25 50 {
		cap drop box
		egen box=group(ESI`j'_gYe ESI`j'_nYe), label
		label list box

		eststo: mlogit box c.age##c.age##c.age##c.age##i.ccode if inrange(age,16,85)
		predict b00_`j' b01_`j' b10_`j' b11_`j' if e(sample), pr
	}

	preserve
		keep if inrange(age,16,85)

		collapse (mean) b??_??, by(ccode age)
		drop if ccode==.
		
		reshape long b00_ b01_ b10_ b11_ , i(ccode age) j(dropsize)
		rename b??_ b??
		label var b00 "MI:0, DI: 0"
		label var b01 "MI:0, DI: 1"
		label var b10 "MI:1, DI: 0"
		label var b11 "MI:1, DI: 1"
		
		gen ESImi_=b10+b11
		label var ESImi_ "25% drop (MI)"
		gen ESIdi_=b01+b11
		label var ESImi_ "25% drop (DI)"
		gen rr_=b10/(b10+b11)
		label var rr "Risk reduction"
		
		kountry ccode, from(cown) to(iso3c)
		rename _ISO3C iso3c

		drop ccode
		rename b?? b??_ 

		reshape wide *_, i(age dropsize) j(iso3c) str
		
		foreach c in AUS CHE DEU GBR USA {
			foreach dv in b00_ ESImi_ ESIdi_ rr_ {
				label var `dv'`c' "`c'"
			}
		}
		
		note: foreach j in 10 25 50 {
		note: egen box=group(ESI`j'_gYe ESI`j'_nYe), label
		note: eststo: mlogit box c.age##c.age##c.age##c.age##i.ccode if inrange(age,16,85)
		note: predict b00_`j' b01_`j' b10_`j' b11_`j' if e(sample), pr
		note: ./resultsets/figure_2_6.do
		note: Created on `= c(current_date)'
		saveold ./resultsets/figure_6.dta, replace

		foreach j in 10 25 50 {
			foreach dv in b00_ ESImi_ ESIdi_ rr_ {
				if "`dv'"=="b00_"   local t "No `j'% drop"
				if "`dv'"=="ESImi_" local t "`j'% drop in MI"
				if "`dv'"=="ESIdi_" local t "`j'% drop in DI"
				if "`dv'"=="rr_"    local t "Risk reduction (`j'% drops)"
				tw (connected `dv'* age if dropsize==`j', sort ms(none none none none none) ), legend(r(1) pos(6)  region(lstyle(none)) ) /// lp(solid dash dot shortdash shortdash_dot )
					xlabel(15(10)85, angle(0)) ///
					ytitle(`t') xtitle("") ///
					scheme(rbn1mono)
				more
				graph save ./resultsets/figures/tmp/tmpFig6_`dv'`j'.gph, replace
			}
		}	
		
		foreach j in 10 25 50 {
			grc1leg ///
				./resultsets/figures/tmp/tmpFig6_ESImi_`j'.gph ///
				./resultsets/figures/tmp/tmpFig6_ESIdi_`j'.gph ///
				./resultsets/figures/tmp/tmpFig6_b00_`j'.gph ///
				./resultsets/figures/tmp/tmpFig6_rr_`j'.gph, scheme(rbn1mono)
			
			graph save   ./resultsets/figures/figure6_`j'.gph, replace
			graph export ./resultsets/figures/figure6_`j'.pdf, replace
			*graph export ./resultsets/figures/figure6_`j'.emf, replace
		}
		
	restore		

display "$S_TIME  $S_DATE"
cap log close

